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Particle shape dependence in 2D granular media 
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Abstract - Particle shape is a key to the space-filling and strength properties of granular matter. 
We consider a shape parameter n describing the degree of distortion from a perfectly spherical 
shape. Encompassing most specific shape characteristics such as elongation, angularity and non- 
convexity, n is a low-order but generic parameter that we used in a numerical benchmark test for a 
systematic investigation of shape-dependence in sheared granular packings composed of particles 
of different shapes. We find that the shear strength is an increasing function of n with nearly the 
same trend for all shapes, the differences appearing thus to be of second order compared to g. We 
also observe a nontrivial behavior of packing fraction which, for all our simulated shapes, increases 
with g from the random close packing fraction for disks, reaches a peak considerably higher than 
that for disks, and subsequently declines as g is further increased. These findings suggest that 
a low-order description of particle shape accounts for the principal trends of packing fraction 
and shear strength. Hence, the effect of second-order shape parameters may be investigated by 
considering different shapes at the same level of g. 



k"* The hard-sphere packing is at the heart of various mod- 
k> els for the rheology and (thcrmo) dynamical properties of 
amorphous states of matter including liquids, glasses and 
Agranular materials [UH]. Such models reflect both the 
purely geometrical properties of sphere packings, e.g. the 
order-disorder transition with finite volume change [3], 
and emergent properties arising from collective particle in- 
teractions, e.g. force chains and arching in static piles [3]. 
As to non-spherical particle packings, rather recent results 
suggest that such packings exhibit higher shear strength 
than sphere packings 5 15 , and may approach unusually 
high packing fractions pi fTBT - HS] . However, a systematic 
and quantitative investigation of shape-dependence is still 
largely elusive since particle shape characteristics such as 
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elongation, angularity, slenderness and nonconvexity are 
described by distinct groups of parameters, and the effect 
of each parameter is not easy to isolate experimentally. 

In order to evaluate the shape-dependence of gen- 
eral granular properties such as packing fraction, shear 
strength and internal structure for particles of different 
shapes, we designed a numerical benchmark test that was 
simulated and analyzed by the members of a collabora- 
tive group (CEGEO). The idea of this test is that various 
non-spherical or non-circular shapes can be characterized 
by their degree of distortion from a perfectly spherical or 
circular shape. Let us consider an arbitrary 2D shape as 
sketched in Fig. [TJ The border of the particle is fully 
enclosed between two concentric circles: a circumscribing 
circle of radius R and an inscribed circle of radius R — AR. 
We define the n-set as the set of all shapes with borders 
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enclosed between a pair of concentric circles (spheres in 
3D), touching both circles and having the same ratio 



Four different particle shapes belonging to the same 77- 
set are shown in Fig. [2J A non-zero value of 77 corresponds 
to non-convexity for A-shape, elongation for B-shape, an- 
gularity for C-shape, and a combination of angularity and 
elongation for D-shape. 




Fig. 2: Four different shapes belonging to the same 77-set with 
77 = 0.4: trimer (A), rounded-cap rectangle (B), truncated tri- 
angle (C), and elongated hexagon (D). 




Fig. 3: Snapshots of the simulated packings in the densest 
isotropic state for 77 = 0.4. 

The parameter 77 is obviously a rough low-order shape 
parameter; see also [T^]. But, encompassing most spe- 
cific shape parameters, it provides a general framework in 



Fig. 4: Shear strength sin tp* of packings composed of various 
particle shapes (see Fig. |2) as a function of 77. 



which shape-dependence may be analyzed among parti- 
cles of very different shapes. Within an 77-set, each spe- 
cific shape may further be characterized by higher-order 
parameters. The issue that we address in this Letter is to 
what extent the packing fraction and shear strength are 
controlled by 77 and in which respects the behavior depends 
on higher-order shape parameters . 

The benchmark test is based on the four shapes of Fig. 
[2] The A-shape (trimer) is composed of three overlapping 
disks touching the circumscribing circle and with their in- 
tersection points lying on the inscribed circle; the B-shape 
(rounded-cap rectangle) is a rectangle touching the in- 
scribed circle and juxtaposed with two half-disks touching 
the circumscribing circle; the C-shape (truncated trian- 
gle) is a hexagon with three sides constrained to touch the 
inscribed circle and all corners on the circumscribing cir- 
cle; and the D-shape (elongated hexagon) is an irregular 
hexagon with two sides constrained to touch the inscribed 
circle and two corners lie on the circumscribing circle. The 
range of geometrically defined values of 77 for a given shape 
(defined by a construction method) has in general a lower 
bound 770. For A and B, the particle shape changes con- 
tinuously from a disk, so that 770 = whereas we have 
7/o = l- \/3/2 ^ 0.13 for C and D. 

Two different discrete element methods (DEM) were 
used for the simulations: contact dynamics (CD) and 
molecular dynamics (MD). In the CD method, the par- 
ticles are treated as perfectly rigid [20] whereas a linear 
spring-dashpot model was used in MD simulations with 
stiff particles (k n /po > 10 3 , where k n is the normal stiff- 
ness and po refers to the confining pressure) [21]. The 
trimers were simulated by both methods for all values of 
77. We refer below as A (for CD) and A' (for MD) to these 
simulations. The packing C was simulated by MD whereas 
the packings B and D were simulated by CD. In CD sim- 
ulations, the coefficient of restitution was set to zero. In 
MD simulations, the damping parameter was taken very 
close to the critical damping coefficient so that the resti- 
tution coefficient was also negligibly small [22]. Note that 
in quasi-static flow, the relaxation time of the particles is 
short enough (compared to the inverse shear rate) to al- 
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Fig. 5: Friction mobilization in the steady state as a function Fig. 6: Packing fraction in the isotropic state as a function of 
of r) for different particle shapes. V for different particle shapes. 



low for efficient dissipation of kinetic energy in each time 
step. For this reason, in contrast to granular gases, the 
exact values of the damping parameters or restitution co- 
efficients have practically no influence on the numerical 
data analyzed below |23j . 

For each shape, several packings of 5000 particles were 
prepared with r\ varying from to 0.5. To avoid long- 
range ordering, a size polydispersity was introduced by 
taking R in the range [R m in, Rmax] with R max = ?>R min 
and a uniform distribution of particle volumes. A dense 
packing composed of disks (77 = 0) was first constructed by 
means of random deposition in a box (24] . For other values 
of ?7, the same packing was used with each disk serving 
as the circumscribing circle. The particle was inscribed 
with the desired value of 77 and random orientation inside 
the disk. This geometrical step was followed by isotropic 
compaction of the packings inside a rectangular frame. 
The gravity g and friction coefficients between particles 
and with the walls were set to during compaction in order 
to avoid force gradients. Fig. [3] displays snapshots of the 
packings for rj = 0.4 at the end of isotropic compactiorQ. 

The isotropic samples were sheared by applying a slow 
downward velocity on the top wall with a constant con- 
fining stress acting on the lateral walls. During shear, the 
friction coefficient /j between particles was set to 0.5 and 
to with the walls. The shear strength is characterized 
by the internal angle of friction ip defined by 
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where the subscripts 1 and 2 refer to the principal stresses. 
shop increases rapidly from zero to a peak value before 
relaxing to a constant material-dependent value srnip*, 
which defines the shear strength at large strain at a steady 
stress state. 

Figure 2] shows the dependence of sin ip* with respect 
to rj for our different shapes. Remarkably, sin tp* increases 
with 77 at the same rate for all shapes. The data nearly co- 
incide between the A and B shapes, on the one hand, and 



between C and D shapes, on the other hand. This sug- 
gests that nonconvex trimers and rounded-cap rectangles, 
in spite of their very different shapes, belong to the same 
family (rounded shapes). In the same way, the truncated 
triangles and elongated hexagons seem to belong to the 
family of angular particles and exhibit a shear strength 
slightly above that of rounded shapes. Note also that the 
results are robust with respect to the numerical approach 
as the packings A and A' were simulated by two different 
methods. 

The increase of shear strength with 77 may be attributed 
to the increasing frustration of particle rotations as the 
shape deviates from a disk [TT][25]. Since the particles 
may interact at two or three contact points (A-shape) or 
through sidc-to-side contacts (shapes B, C and D), the 
kinematic constraints increase with 77 and frustrate the 
particle displacements by rolling. The restriction of rolling 
leads to enhanced role of friction in the mechanical equilib- 
rium and relative sliding of particles during deformation. 
A related static quantity is the mean friction mobilization 
defined by M = (ft/(lJ-fn)}, where f t is the magnitude of 
the friction force, /„ is the normal force, and the average 
is taken over all force-bearing contacts in the system. 

To evaluate the effect of particle shape, we consider the 
parameter 

M(t?) 



M„ 



M(t] = 0) 



1 



(3) 



1 Animation videos of the simulations can be found at www.cgp- 
gat c way. org/ref012. 



as a function of 77 for different shapes, where M(r] = 0) 
is the friction mobilization for circular particles. Fig. [5] 
shows that M v is a globally increasing function of 77 for all 
shapes. The parameter 77 appears also in this respect to 
account for the global trend of friction mobilization, and 
the differences observed in Fig. [5] among different shapes 
are rather of second order. 

We also observe that the proportions of double and 
triple contacts for A-shapc packings and the proportion of 
side-to-side contacts for other shapes increase with 77. For 
noncircular particles, one should distinguish the coordina- 
tion number Z, defined as the mean number of contact- 
ing neighbors per particle, from the "contact coordination 
number" Z c defined as the mean number of contacts per 
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Fig. 7: Pore volume reduction by (a) overlap between self- 
porosities; (b) steric pores. 



particle. Obviously, for the calculation of both Z and Z c 
only the force-bearing contacts and non-floating particles 
are taken into account [55]. We have Z = Z c ~ 4 for 
the disks in the initial state prepared with p — 0. This 
value corresponds to an isostatic state in which one ex- 
pects Z = 2Nf, where Nf is the number of degrees of 
freedom of a particle [57] ■ For frictionless disks, we have 
Nf = 2 (two translations] degrees of freedom) , leading to 
Z = 4. For noncircular shapes, we have Nf = 3 since the 
rotational degrees of freedom take part in the mechanical 
equilibrium of the particles. Hence, if isostaticity holds 
also for noncircular frictionless particles, we expect Z = 6. 
We observe instead Z < 5 for all our packings. However, 
we find Z c ~ 6 for r\ 7^ if each side-to-side contact is 
counted twice, representing two independent constraints. 
This result is consistent with the isostatic nature of a pack- 
ing of frictionless noncircular particles and shows that the 
packings of noncircular shape are not under-constrained as 
previously suggested [55]. For p = 0.5, the packings are 
no more isostatic and Z and Z c vary only slightly with r\ 
with values in the range 3 to 4 for Z and in the range 4 
to 5 for Z c in the course of shearing. 

We now focus on the packing fraction which crucially 
depends on particle shape. Fig. [5] shows the packing frac- 
tion p %s ° in the initial isotropic state as a function of 77. 
We observe a nontrivial behavior for all particle shapes: 
the packing fraction increases with 77, passes by a peak de- 
pending on each specific shape and subsequently declines. 
For the B-shape a sharp decrease of p lso occurs beyond 
77 = 0.5 as was shown in |10j . 

This unmonotonic behavior of packing fraction was 
observed by experiments and numerical simulations for 
spheroids as a function of their aspect ratio [2i rT6lfT7lf28l - 
150] . The decrease of the packing fraction is attributed to 
the excluded-volume effect that prevails at large aspect 
ratios and leads to increasingly larger pores which cannot 
be filled by the particles [5S] . The observation of this un- 
monotonic behavior as a function of 77 for different shapes 
indicates that it is a generic property depending only on 
deviation from circular shape. This behavior may thus be 
explained from general considerations involving the pa- 
rameter 77 but with variations depending on second-order 
shape characteristics. 
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Fig. 8: Normalized packing fractions fitted by Eq. (0. 



A plausible second-order parameter is 
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nR 2 ' 
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where V p is the particle volume in 2D. Its complement 
I — v is the "self-porosity" of a particle, i.e. the unfilled 
volume fraction inside the circumscribing circle. Keeping 
the radius R of the circumscribing circle constant, p lso = 
Vp/V varies with 77 as a result of the relative changes of 
V p and the mean volume V per particle. The free (pore) 
volume per particle is Vf = V — V p . 

At 77 = 0, the free volume Vf is only composed of steric 
voids, i.e. voids between three or more particles, and the 
packing fraction is given by p(0) = ttR 2 /V(0). For 77 > 0, 
the void patterns are more complex but can be described 
by considering the generic shape of particles belonging to 
a given 77-set. The borders of a particle involve "hills", 
which are the parts touching the circumscribing circle, and 
"valleys" touching the inscribed circle. The volume V per 
particle varies with 77 by two mechanisms. First, the hills 
of a particle may partially fill the valleys of a neighboring 
particle; Fig. [7Ja). Secondly, the steric voids between the 
hills shrink as 77 increases due to the increasing local curva- 
ture of the touching particles; Fig. UJb). To represent this 
excess or loss of pore volume due to the specific jamming 
configurations induced by particle shapes, we introduce 
the function h(rj) by setting 



v(v) = Vim) - 7ri? 2 /i(7/0, 
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with h(r)o) = 0. With these assumptions, the packing 
fraction is expressed as 



1 - h(r])p{r) )' 



(6) 



The function v{rf) is known for each shape but h{rj) 
needs to be estimated. A second-order polynomial ap- 
proximation 



Hv) = a (v - Wo) + P{w - Wo? 



(7) 



together with Eq. [5] allows us to recover the correct trend 
and to fit the data as shown in Fig. \E\ The error bars 
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represent the variability at 770 assumed to be the same for 
all other values of rj. The parameter a ensures the increase 
of packing fraction with 77 at low values of the latter and it 
basically reflects the shrinkage of steric pores (Fig. E£b)) 
whereas /3 accounts for the overlap between circumscribing 
circles (Fig. EJa))) and is responsible for the subsequent 
decrease of the packing fraction. 

The fitting parameters in Fig. \8\ arc a ~ 1.30, 1.29, 
1.14, 1.17 and j3 ~ 1.23, 1.20, 0.23, 0.20 for C, A, D and B 
shapes, respectively with increasing peak value. Note that 
the values of /3 are considerably smaller for B and D that 
have an elongated aspect and for which the overlapping of 
self-porosities prevails as compared to A and C for which 
the shrinkage of the initial pores is more important. 

In summary, our benchmark simulations show that a 
low-order shape parameter 77, describing deviation with 
respect to circular shape, controls to a large extent both 
the shear strength and packing fraction of granular me- 
dia composed of noncircular particles in 2D. The shear 
strength is roughly linear in 77 whereas the packing frac- 
tion is unmonotonic. Our simple model for this unmono- 
tonic behavior is consistent with the numerical data for 
all shapes. It is governed by a first-order term in 77 for 
the shrinkage of the initial steric pores and a second-order 
term in 77 for the creation of large pores by shape-induced 
steric pores. The effect of higher-order shape parameters 
may be analyzed also in this framework in terms of differ- 
ences in packing fraction and shear strength among various 
shapes belonging to the same 77-set. An interesting issue 
to be addressed in future is whether a generic second-order 
parameter accounting for such differences exists. Another 
aspect that merits further investigation is the joint ef- 
fects of size polydispersity and particle shape. The shear 
strength is independent of particle size polydispersity as a 
result of the capture of force chains by the class of larger 
particles [3T] . But the packing fraction and force and con- 
tact anisotropy depend on both shape and polydispersity. 
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